Structure and relaxations in liquid and amorphous Selenium 
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We report a molecular dynamics simulation of selenium, described by a three-body interaction. 
The temperatures T g and T c and the structural properties are in agreement with experiment. The 
mean nearest neighbor coordination number is 2.1. A small pre-peak at about 1 A -1 can be explained 
in terms of void correlations. In the intermediate self-scattering function, i.e. the density fluctuation 
correlation, classical behavior, a- and /3-regimes, is found. We also observe the plateau in the (3- 
regime below T s . In a second step, we investigated the heterogeneous and/or homogeneous behavior 
of the relaxations. At both short and long times the relaxations are homogeneous (or weakly 
heterogeneous). In the intermediate time scale, lowering the temperature increases the heterogeneity. 
We connect these different domains to the vibrational (ballistic), j3- and a-regimes. We have also 
shown that the increase in heterogeneity can be understood in terms of relaxations. 



PACS numbers: 61.20.Lc, 66.10-x, 61.43. Dq, 61.43.Fs 



I. INTRODUCTION 

Although glass is one of the most common materials, 
the glass transition is still poorly understood. It is a 
continuous transition in which the viscosity of the glass 
forming liquid increases from 10 -3 Pa s in the liquid to 
more than 10 9 Pa s in the super-cooled state. It is, there- 
fore, easy to understand that very different time scales 
become important near the glass transition, and different 
types of relaxations might be observed.- 

The Mode Coupling Theory (MCT)W gives a micro- 
scopic picture of this transition. This theory focuses on 
the density correlation function &(q, t), the intermedi- 
ate self-scattering function, and proposes a mechanism of 
back- flow to explain the increase of the viscosity and/or 
of the time-scalescl. One of the most striking results of the 
MCT is the prediction of a critical temperature T c below 
which the system becomes non-ergodic. In other words 
the system is trapped in a well of the energy landscape. 
This feature is related to a non-zero value of $((?, t) when 
t — > co. Above T c , the function <j>(g, t) shows a short time 
relaxation, related to the vibrational (often called ballis- 
tic) regime, and a long time one, also called ev-relaxation. 
Below T c , a third regime appears, the so-called /3-regime, 
which is first seen as a shoulder and saturates at a finite 
value below T„. 

This non-ergodieity has been seen in many experimen- 
tal measurementstaEl and computer simulationsuiI9 on 
different types of materials ranging from the fragile poly- 
meric glasses to strong glasses, such as Si02- In this 
paper we want to go a step further. Using a model of 
selenium, we check for the presence of these two or three 
(depending on the temperature) relaxation steps, and ask 
the following question. Does each atom have the same 
probability to relax? If below T c the system becomes 
non-ergodic, and is trapped in a well of the energy land- 
scape, are all atoms still equivalent, or are some more (or 
less) mobile than others? We can reformulate this ques- 
tion and ask whether the relaxations are homogeneous 



or heterogeneous. According to some authorsOEJ the 
relaxations should be more horoogpjeous, in particular 
in the cv-regime, whereas othersliiaEil claim that the re- 
laxations in amorphous or disordered materials are more 
heterogeneous. The answer therefore does not seem to 
be obvious. 

From the theoretical point of view, simple one-atomic 
systems such as soft or hard spheres or Lennard- Jones 
systems would be optimal to study. Unfortunately these 
simple systems crystallize rapidly at temperatures near 
the glass transition temperature (T g ) and, therefore, can 
be utilized only for studies in the liquid, well above T g , or 
deep in the glassy state, T <C T g . One possibility to avoid 
crystallization is to introduce special features in the inip;- 
atomic interaction potential which penalize ordcringEat 3 ]. 
The most common approach is to use binary mixtures of 
atoms. A different approach is to simulate a real mono- 
atomic glass former, such as selenium, which has , 
studied extensively in experiment, see the reviewsc 
Se is covalently bound and prefers a coordination num- 
ber of two. This is reflected in the different crystal struc- 
tures. The most stable trigonal phase consists of parallel 
helical chains. Two monoclinic forms are composed of 
rings of eight atoms. The polymorphs are distinguished 
by the correlation between neighboring dihedral angles. 
Depending on this correlation one has either a trans- 
chains) or a cis-configuration (rings). The energy dif- 
ference between the cis- and trans-configuration was es- 
timated as only 0.03 eVcJ. This low energy difference 
implies that in glasses both configurations should coex- 
ist, which in turn strongly hinders crystallization. From 
a first-principles moleculax-. dynamics simulation using 
64 atoms Hohl and JonescJ conclude that both amor- 
phous and liquid selenium can be viewed as consisting of 
branched chains which include rings of different length. 
The small size prevented quantitative statistics of chain 
and ring lengths and branching points. The fraction of 
atoms having twofold cpoirdjnation varies in the literature 
between 95% and 71%ErH. 
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To study dynamical properties, larger systems are 
needed and one has to resort to effective inter-atomic 
interactions. This immediately leads to the difficulty of 
simultaneously having to describe the covalent binding 
in the chains and rings and the van-der-Waals like inter- 
action between the rings, as well as possible branching 
and bond breaking. One possibility is to disregard the 
last two, and to use different interactions for atoms in 
the same chain and in different chains, respectively. Sim- 
ilar to simulations of polymers one then considers a glass 
or a liquid of chains of a predefined length. This frag- 
mented chain method has been employed both for elec- 
tronic structure calcuiatianscd and for classical molecular 
dynamics simulationsE3'E3. Another possibility is to use a 
more simple description-af the electronic properties, such 
as tight binding modekc3. 

We follow a different line. As done before by Stillinger 
et al. in their study of liquid sulfma, we use one effective 
inter-atomic potential for both, the intra-chain and the 
inter-chain interactions. 

This paper is laid out as follows: in section II we re- 
port the details of the simulations, and of the production 
of the liquid and glassy samples used in the measure- 
ments of the different properties reported in this work. 
Section III is devoted to the determination of the glass 
transition temperature T g and the critical temperature 
T c . Given these temperatures, we report the evolution of 
the structure through the glass transition in section IV. 
In the next section we focus on the relaxations and the 
intermediate self-scattering function. In section VI we 
present the tools, used to study the homogeneity or het- 
erogeneity of these relaxations, report the measurements 
and discuss them. Finally we conclude. 



II. SIMULATIONS 

We performed classical molecular dynamics simula- 
tions on a system of N = 2000 atoms interacting via a 3- 
body potential. This potential was built to reproduce the 
properties of small clusters of selenium and of the trigo- 
nal crystalline phase. Details of the potential are given in 
Ref.E3. The potential has previously been used to calcu- 
late,the vibrationsEj and local relaxations in amorphous 
SeE3. In these simulations it was shown that both the low 
frequency resonant vibrations and the local relaxations 
are centered on groups of ten and more atoms. This 
finding is in agreement with the interpretation of experi- 
ments by the soft potential modelES. From a Monte Carlo 
study of liquid Se it was concluded that the model pro- 
vides a sound basis for the study of both the microscopic 
and the electronic structure, despite some deficiency in 
the treatment of the van-der-Waals interactional 

During the simulations the pressure was fixed to zero 
Pa, i.e. we worked with equilibrium structures. In or- 
der to keep the prpsswe constant we used a Parrincllo- 
Rahman algorithme3H. The temperature was kept con- 



stant by rescaling the velocities at each integration step. 
We checked that the way we control the pressure and 
temperature influenced neither the dynamics of the sys- 
tem nor the results of our simulations. 

The equations of moiipn were integrated using the ve- 
locity Verlet algorithrro. Taking care of the stability of 
the algorithm, we chose the time steps equal to 1 fs in 
the liquid, 2 fs in the glassy phase and 4 fs for the lowest 
temperature (6 K). 

To improve the statistics we used four independent 
starting configurations to obtain the different samples 
used in the measurement. These samples were produced 
as follows: first we equilibrated a liquid at 550 K (above 
the melting point) , and cooled it to 290 K with a quench 
rate of 10 13 K/s. At this temperature we let the systems 
equilibrate during 8 ns and then quenched them to K 
with the same quench rate. During both quenches we 
saved configurations at several temperatures and again 
equilibrated them before using them as input of the mea- 
surements. The equilibration times were 8 ns above 
290 K, 16 ns between 290 K and 6 K, and 32 ns at 6 K. 
After the equilibration period several relevant physical 
quantities were computed. 



III. 



DETERMINATION 



To obtain the relevant temperature scale, we first de- 
termined the glass transition temperature T g . For this we 
followed the evolution of the volume during the quench 
process. As the coefficient of volume expansion is higher 
in the liquid than in the solid phase, one observes a 
change of slope of the volume curve when the system 
is quenched through the glass transition. 




FIG. 1. Evolution of the atomic volume of liquid and amor- 
phous Se atoms during the quench. The change of slope be- 
tween high and low temperatures determines the glass tran- 
sition temperature T g . 

From Fig. 1 the glass transition temperature is es- 
timated, as T g » 300 K. Experimentally T g is about 
305 KEI. The good agreement between the simulated 
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and experimental values of T g might seem surprising re- 
garding the high quench rate and the usual discrepancies 
observed in simulation. However, one should notice that 
due to the aging over several ns the effective quench rate 
is lower, » 10 10 K/s. 

Another characteristic temperature is the critical tem- 
perature T c given by the Mode Coupling Theory. This 
temperature can be obtained from the diffusion constant 
D which- according to the MCT, follows a power law 
(T— T c ) 7 El. The diffusion constant is related to the atomic 
mean square displacement by 



D = lim 



(|r(i + M-r(* )| 2 > 



(if 



(1) 



The diffusion constants obtained, Fig. 2, are in .rcason- 
able agreement with experiments on liquid SeCJCJ. At 
the higher temperatures D can be fitted by an Arrhc- 
nius law with an activation energy of 0.3 eV, in agree- 
ment with results from first-principle molecular dynam- 
ics simulationsE3. In the inset of Fig. 2 we also show, 
by a dashed line, the fit by the MCT power law. Due 
to the correlation between T c and 7 such a fit can only 
approximately determine these values. 





* 




10' 








o'° ^ - 

c- 

,-">' 


10' 

10" f 

10' 

f 

10- 


10' 10 2 

T-T . (K) 


/ 

f 

1 




i 

4 




/ 

y 



100 200 300 400 500 600 

T(K) 

FIG. 2. Diffusion constant of liquid and amorphous Se as 
function of temperature. The full circles show the values ob- 
tained from the MD simulations. The dashed line shows a fit 
with the power law D oc (T — Tc) 1 ■ The inset shows the dif- 
fusion constant D (o) and the decay time of the a-relaxation 

in a log-log representation, the y-axis 
-l 



t- 1 (o) versus (T-T c 
is in units of 10~ 10 m 2 



1 for D, and ps 1 for r 



Fixing T c about 10% above T s , i.e. T c = 330 K we 
get 7 = 1.88. From the same fit to the decay time of 
the a-process, see section V, we obtain for T c = 330 K a 
value 7 = 1.86, which is in excellent agreement. 



IV. STRUCTURAL PROPERTIES 

Having obtained the relevant temperature scales we 
now turn to structural properties. First we compute the 
pair correlation function (PCF) at temperatures ranging 



from the liquid down into the glass. The PCF is defined 
by: 




9(r) 



where (...) denotes the average over configurations. 



(2) 
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FIG. 3. Pair correlation function of Se at three different 
temperatures: 6 K (solid line), 290 K (dotted line), 550 K 
(long dashed line). The inset shows the mean coordination 
numbers for the same temperatures. 

Fig. 3 shows the PCF for three different temperatures, 
550 K in the liquid phase, 290 K just below T g and 6 K 
deep in the glassy phase. The inset shows the mean co- 
ordination number. The positions of the peaks do not 
change strongly upon cooling, but broaden markedly. As 
usual oscillations at large distances are more strongly 
damped at high temperatures. The spatial correlations 
at large distances weaken with increasing temperature. 
In all cases the correlations are weak for distances greater 
than 10 A. The minimum near 4 A for low temperatures 
signals that the second neighbor shell becomes more pro- 
nounced. The mean coordination number (Fig. 3 inset) 
indicates the same behavior. The mean nearest neighbor 
coordination is about at all temperatures, similar to 
the experimental valucc3. This value of around 2 is the 
signature of the chains and rings forming the amorphous 
selenium structure. The small excess of 0.1 compared 
to the ideal value of 2 indicates a prevalence of over- 
coordinated atoms (branching) over under-coordinated 
ones (chain ends). At the lowest temperature (6 K) we 
also observe a small plateau in the coordination number 
near 4 A. The change of neighbor numbers with tempera- 
ture for larger distances reflects the lower density at high 
temperatures. This indicates that with increasing tem- 
perature the chain structure remains, but the distance 
between chains increases. 

From the PCF the structure factor S(q) can be com- 
puted by a spatial Fourier transform: 
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S(q) = 1 + 1 r^{g{r) l)^dr (3) 
N J qr 



Fig. 4 shows S(q) for the three temperatures used 
in Fig. 3. As in the PCF, the peaks do not shift 
strongly with temperature, they merely become more 
damped with increasing temperature. The positions of 
the peaks agr ee . wi th experimentsoEj and previous sim- 
ulations on SeE3ca. In addition we see a small pre-peak 
near q — 1 A -1 . Experimentally the pre-peak in amor- 
phous selenium merges with the first diffraction peak at 
about 2 A -1 and is only seen as a shoulder. 
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FIG. 4. Structure factor of Se at the same temperatures 
as Fig. 3: 6 K (solid line), 290 K (dotted line), 550 K (long 
dashed line). 

To study this pre-peak, we quenched two more sets 
of 10 independent liquids, each to zero K applying two 
different pressures: zero pressure and 10 GPa. Finally 
we minimized the potential energy for both sets, and 
released the pressure for the second set. This gave us 
at T = K two sets of samples with different densities 
p = 4400 kg/m 3 and p — 5090 kg/m 3 , both with equilib- 
rium structures. The average potential energy per atom 
of the low density samples is only 3.5 meV less than the 
one at the high density. This low value might indicate 
that at high temperature voids are present in thermody- 
namic equilibrium. For both sets of equilibrium struc- 
tures we computed the structure factors by 



S(g) = /^exp(iq(r i (t)-r i (*)))\ , (4) 



where q are the q- vectors compatible with the simulation 
box. The minimal g-values, given by the periodic bound- 
ary conditions are q = 0.17 A" 1 and q = 0.16 A -1 for 
the high and the low density samples, respectively. 
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FIG. 5. Structure factor of Se at T = K for two different 
densities: p = 4400 kg/cm 3 (solid line) and p — 5090 kg/cm 3 
(dashed line). 

Whereas the low density samples show a small pre- 
peak near lA _1 , it is absent in those of high density 
(Fig. 5). Checking the mean coordination number at 
the two densities (Fig. 6), one sees that the number of 
first neighbors changes very little with density: there are 
chains and rings at both densities. The number of second 
neighbors, however, is reduced for the lower density. 
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FIG. 6. Mean coordination number of Se at T — K at 
two different densities: p — 4400 kg/cm 3 (solid line) and 
p = 5090 kg/cm 3 (dashed line). 

This is the same effect as previously noted for the tem- 
perature dependence. When the density is low, i.e. when 
the volume is high the system is constituted of chains and 
rings which are further apart from each other. In other 
words we interpret the pre-peak as the signature of cor- 
relations between holes in the steucture. A similar effect 
was seen in simulations of SiC^t^. As a consequence of 
the too high density of the crystalline structured, the 
density of our amorphous selenium is too high in com- 
parison with the experimental value. 
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V. INTERMEDIATE SELF-SCATTERING 
FUNCTION 

We now focus on the local relaxations. First, we com- 
pute the intermediate self-scattering function (ISSF), the 
correlation function of the density fluctuations, 

$(?,*) = (Sp- q (t + to)Sp q (t )) to . (5) 

This can be rewritten as the spatial Fourier transform of 
the van Hove self correlation function G s (r, t) 

*(«,*) = [°° Gs (r,tf^-dr (6) 
Jo V 

where G s (r, t) is given byS 

G s (r,t) = (8(r-\r i (t + t )-r i (t )\)} i<to . (7) 

The ISSF of Se is not easily accessible in experiment. 
It corresponds to the time Fourier transform of the inco- 
herent dynamic structure factor, but neutron scattering 
by Se is mainly coherent. Nevertheless this quantity is 
accessible to simulation. 




t(ps) 



FIG. 7. Intermediate self-scattering function (ISSF) of Se 
versus time for different temperatures above and below T g . 
From top to bottom: W5K, 200K, 255K, 290 K, 330 K, 
355 K, 400 K, 445 K, 495 K and 550 K. 

In Fig. 7 we report the ISSF at q = 2.1 A -1 corre- 
sponding to the first diffraction peak of the S(q) in Fig. 4. 
A fast decrease of $>(q,t) on the time-scale of a picosec- 
ond, is observed at all temperatures. Decreasing the tem- 
perature a shoulder appears for intermediate time scales. 
As the system approaches the glass transition tempera- 
ture, <E>(g, t) starts to saturate and finally shows a plateau 
for intermediate and long times. As customary this curve 
is described as follows: First there is the ballistic or vi- 
brational regime (corresponding to the fast decreases at 
low time). Then, for T > T c , $(q,t) goes to zero (the 
so-called a-regime). The shoulder or plateau at lower T 
corresponds to the so-called /3-regime. This plateau in- 
dicates that the system is trapped in a limited area of 
phase-space. 



According to the MCT, the a-regime above T c is driven 
by a master curve which is obtained by rescaling the time 
by a characteristic decay time, r defined as the time when 
the ISSF has drops to 1/e, r) = 1/e. Above T c , simi- 
larly to the diffusion constant D these values r(T) should 
obey a power law r(T) = (T- T c )—*. Fixing T c = 330 K 
we get a good fit with 7 = 1.86, see inset in Fig. 2. Using 
this r the master curve can be written in scaled form asH: 

$(g,Vr)=$ exp(-(Vr)' 3 ). (8) 

Fig. 8 presents the curves for temperatures above T g 
rescaled by their respective r. We get a value (3 = 0.53 
for temperatures ranging from T = 330 K to T = 430 K. 
We do not intend to give an elaborate test of the MCT 
but show the rescaled curves merely to indicate that the 
rescaling procedure holds approximately. 




FIG. 8. Intermediate self-scattering function of Se ver- 
sus the dimensionless variable t/r, where r is defined by 
$(g, t) = 1/e, see inset of Fig. 2. Temperatures from left 
to right: 290 K, 330 K, 355 K, 400 K, 445 K, 495 K and 
550 K. 

The most striking effect is the plateau corresponding 
to the /3-regime. It indicates that the system falls out of 
equilibrium, and that atoms are trapped in a well of the 
energy landscape. This poses the question of whether all 
the atoms are affected equally or not. 

The same question can be asked for the a-regime which 
can be represented by a stretched exponential decay law. 
Such a law can either mean that a stretched exponential 
decay law governs the whole system, or it originates from 
an average of normal exponential laws with different time 
constants. The first picture is called homogeneous sce- 
nario (the system is everywhere the same) and the second 
one heterogeneous. 

VI. HETEROGENEITY OR HOMOGENEITY ? 

To answer this question, we use again the van Hove 
correlation function (VHF) which represents the proba- 
bility for an atom to move a distance r during a time 
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t. If all the atoms have the same mobility the VHF is a 
Gaussian multiplied by the geometrical factor Airr 2 . In 
the following this factor is always thought to be included 
when we speak of Gaussianity. In the opposite case if 
some atoms are trapped or some are more mobile than 
the majority the VHF will, in general, be non-Gaussian. 
As example we show in Fig. 9 the VHF for two differ- 
ent temperatures for the same length range but for two 
different times. Obviously at high temperatures the sys- 
tem has a higher diffusion constant and the atoms will 
move faster over a given distance. More striking is that 
at high temperatures (in the liquid) the VHF has Gaus- 
sian shape, whereas at low temperatures (in the glass) 
the VHF has an extended tail and cannot be described 
by a Gaussian. Some atoms move over much larger dis- 
tances than the average atom, i.e. they have a much 
higher mobility. 




2.0 3.0 
r(A) 

FIG. 9. Van Hove correlation function of Se at two differ- 
ent temperatures (solid line: 550 K; dashed line 255 K) and 
times, t = 2.1 ps for 550 K and t = 36 ns for 255 K, re- 
spectively. The different times reflect the higher mobility (or 
diffusion) in the liquid. 

In order to quantify these findings, and in accordance 
with previous work, ^ue introduce the non-Gaussianity 
parameter (NGP) c^Eil 



a 2 {t) 



3^4 



1, 



(9) 



where [ii and fi^ are the second and fourth moments of 
the VHF, [i2 = ( r2 (t)) an d ^4 = (r 4 (i)), respectively. 
The NGP is identical to zero for a Gaussian VHF. 

Fig. 10 shows, in a log-linear representation, the a 2 
at different temperatures from the liquid to the glass for 
a time span covering 6 decades. At very short times 
the NGP is nearly zero at all temperatures. The limit- 
ing behavior for large times is more difficult to observe, 
especially at low temperatures. Nevertheless the NGP 
clearly tends to zero. Furthermore, all the curves, below 
and above the glass transition, scale to the same curve 
in the short time range as already shown by Kob et al. 
in their study of a binary super-cooled Lennard-Jones 



liquid0. Our work shows that this property persists in 
the glassy phase and for a completely different structure. 
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FIG. 10. Log-linear plot of the non-Gaussianity parameter, 
ct2, of Se versus time for several temperatures. From top to 
bottom : 200 K, 255 K, 290 K, 330 K, 355 K, 400 K, 445 K 
and 495 K. 

In the intermediate time range the NGP has a maxi- 
mum indicating heterogeneity. At high temperatures, in 
the liquid above T c , this maximum is small and located 
around 10 ps. Upon cooling, it slowly moves to higher 
times. For the temperatures below T c , the behavior is dif- 
ferent. The value of the maximum is larger than 2 and 
it grows by a factor 2 upon cooling by 50 K. Similarly 
the position of the maximum shifts by about an order 
of magnitude for each 50 K. These two observations sug- 
gest that as the system is cooled especially below T g , the 
non-Gaussianity becomes more and more pronounced at 
intermediate time scales. 




t(ps) 

FIG. 11. Log-linear plot of the non-Gaussianity parameter 
of Se versus time for two samples (A and B) at low temper- 
ature: 6 K. The insets show the same quantities in a lin- 
ear-linear plot. 



We now focus on the short time behavior at very low 
temperatures. In Fig. 11, we present the evolution of the 



G 



NGP for two different samples A and B (out of our 4 dif- 
ferent samples) at a very low temperature, T = 6 K. The 
inset gives the same curves in linear-linear representation 
to show them clearly at intermediate times. The curves 
coincide during the first ps in the vibrational regime. For 
the larger, intermediate time scale the NGP of sample A 
(Fig. 11. a) oscillates around a value of 0.2, but the one of 
sample B (Fig. ll.b) rises. The two other samples behave 
similarly to sample A. What is the reason of this differ- 
ence in the non-Gaussian behavior of these two kinds of 
samples? The evolution of both total energy and volume 
were equivalent. 
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FIG. 12. Mean square displacement for samples A (top) 
and B (bottom) during the simulation. 

The mean square displacements, however, evolve dif- 
ferently (fig. 12). In sample A (Fig. 12. a) it oscillates 
around a mean value during the entire simulation run, 
whereas it shows steps for sample B (Fig. 12. b). Thus, 
while sample A stays in one well of the energy landscape, 
sample B moves from one well into another. We can iden- 
tify at least two different wells for sample B. We conclude 
that the value «2 ~ 0.2 of the NGP corresponds to the 
vibrations in the liquid and amorphous selenium. Re- 
laxations from one minimum of the energy landscape to 
another lead to an increase in the NGP. It has been shown 
previously that that these 1og*1 relaxations are collective 
jumps of 10 and more atoms.EH The jump length of a sin- 
gle atom is much less than the nearest neighbor distance. 
The different behavior of the samples, shown in Fig. 12 
reflects the low probability for relaxations at low temper- 
ature. It is not restricted to Se but is typical for glasses. 
The same has been reported also from a simulation of 
simple soft sphere glassE^I and is observed experimentally 
in the telegraph noise of the electric resistivity of point 
contactsEl 



VII. CONCLUSION 

In this article we have presented results of a molecu- 
lar dynamics simulation on the structure and relaxations 



of liquid and amorphous Se. We determined the glass 
transition and critical temperatures, the pair correlation 
function and the structure factor. From the pair corre- 
lation function, and in agreement with experiments, we 
concluded that both liquid and amorphous selenium are 
constituted of chains and rings with a mean coordination 
number of 2.1, slightly above the ideal value 2. Rings and 
chains are interconnected. The structure factor shows a 
small pre-peak around 1 A -1 which in experiment only 
shows as a shoulder of the main peak. This pre-peak is 
explained in terms of a correlation of voids between the 
selenium chains. To prove this assumption we computed 
the structure factors of two sets of samples with two dif- 
ferent densities. At the higher density no pre-peak is 
observed. 

The van Hove correlation function was calculated and 
utilized to compute the intermediate scattering function 
and the non-Gaussianity parameter. For the intermedi- 
ate self-scattering function, the time correlation of the 
density fluctuation, we find the classical behavior: at 
short times a rapid decrease corresponding to the bal- 
listic (or vibrational) regime, and at long times a slow 
decay corresponding to the a-regime. When the system 
reaches T c a shoulder and below T g a plateau evolves 
between these two regimes. This corresponds to the (3- 
regime, and to a memory effect of the correlation func- 
tion, in other words the system falls out of equilibrium. 

The non-Gaussianity parameter shows that, both at 
short and long times, the relaxations are homogeneous or 
only weakly inhomogeneous and all the atoms are largely 
equivalent. For the intermediate time range (correspond- 
ing to the /3-regime) NGP depends on the temperature. 
The lower the temperature the higher the NGP, i.e. the 
higher the heterogeneity. At low temperatures the in- 
crease of non-Gaussianity is due to relaxations. We con- 
clude the following scenario for the heterogeneity: at all 
temperatures, both above and below T g , there is a small 
increase of heterogeneity (a2 = 0.2) due to vibrations 
at short time, at intermediate times a pronounced in- 
crease, due to the relaxations especially at temperatures 
below T g , and finally a decrease, due to long range dif- 
fusion (flow motion) . These different regimes correspond 
to the different regimes observed in the intermediate self- 
scattering function &(q,t). This scenario seems to be 
common to different materials. 
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